﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace Solver
{
  public class Solver
  {
    // сетка без выреза.
      public void compute(double[,] v, double[,] f, double h, double k, int n, int m, out int numStep, int maxStep, out double currentEps, double epsMin)
      {
        double vOld, vNew;
        double h2 = 1.0 / (h * h);
        double k2 = 1.0 / (k * k);
        double a2 = -2 * (h2 + k2);

        currentEps = 0;
        numStep = 1;

        double[,] r = new double[n + 1, m + 1];

        double[,] hs = new double[n + 1, m + 1];


        do
        {
          double curEps = currentEps = 0;




          for (int j = 1; j < m; ++j)
          {
            for (int i = 1; i < n; ++i)
            {
              r[i, j] = (f[i, j] + (v[i + 1, j] - 2 * v[i, j] + v[i - 1, j]) * h2 + (v[i, j + 1] - 2 * v[i, j] + v[i, j - 1]) * k2);

            }

          }

          if (numStep == 1)
          {
            for (int j = 1; j < m; ++j)
            {
              for (int i = 1; i < n; ++i)
              {
                hs[i, j] = -r[i, j];
              }
            }
          }
          else
          {
            double betta1 = 0;
            double betta2 = 0;
            for (int j = 1; j < m; ++j)
            {
              for (int i = 1; i < n; ++i)
              {
                betta1 += ((hs[i + 1, j] - 2 * hs[i, j] + hs[i - 1, j]) * h2 + (hs[i, j + 1] - 2 * hs[i, j] + hs[i, j - 1]) * k2) * r[i, j];
                betta2 += ((hs[i + 1, j] - 2 * hs[i, j] + hs[i - 1, j]) * h2 + (hs[i, j + 1] - 2 * hs[i, j] + hs[i, j - 1]) * k2) * hs[i, j];
              }
            }
            double betta = betta1 / betta2;

            for (int j = 0; j < m + 1; ++j)
            {
              for (int i = 0; i < n + 1; ++i)
              {
                hs[i, j] = -r[i, j] + betta * hs[i, j];
              }
            }
          }

          double alpha1 = 0;
          double alpha2 = 0;
          for (int j = 1; j < m; ++j)
          {
            for (int i = 1; i < n; ++i)
            {
              alpha1 += r[i, j] * hs[i, j];
              alpha2 += ((hs[i + 1, j] - 2 * hs[i, j] + hs[i - 1, j]) * h2 + (hs[i, j + 1] - 2 * hs[i, j] + hs[i, j - 1]) * k2) * hs[i, j];
            }
          }

          double alpha = -alpha1 / alpha2;
          if (Math.Abs(alpha) > 1e-12)
          {
            for (int j = 1; j < m; ++j)
            {
              for (int i = 1; i < n; ++i)
              {
                vOld = v[i, j];

                vNew = v[i, j] + alpha * hs[i, j];

                curEps = Math.Abs(vOld - vNew);
                if (curEps > currentEps)
                {
                  currentEps = curEps;
                }
                v[i, j] = vNew;

              }
            }
            numStep++;
          }

        } while (currentEps > epsMin && numStep <= maxStep);
        numStep--;
    }

    // сетка с вырезом.
    // вырез (1/2, 1/2) левый нижний
      public void compute2(double[,] v, double[,] f, double h, double k, int n, int m, out int numStep, int maxStep, out double currentEps, double epsMin)
      {
        double vOld, vNew;
        double h2 = 1.0 / (h * h);
        double k2 = 1.0 / (k * k);
        double a2 = -2 * (h2 + k2);

        currentEps = 0;
        numStep = 1;

        double[,] r = new double[n + 1, m + 1];

        double[,] hs = new double[n + 1, m + 1];


        do
        {
          double curEps = currentEps = 0;




          for (int j = 1; j < m; ++j)
          {
            for (int i = 1; i < n; ++i)
            {
              if ((j <= m / 2) && (i <= n / 2))
                continue;

              r[i, j] = (f[i, j] + (v[i + 1, j] - 2 * v[i, j] + v[i - 1, j]) * h2 + (v[i, j + 1] - 2 * v[i, j] + v[i, j - 1]) * k2);

            }

          }

          if (numStep == 1)
          {
            for (int j = 1; j < m; ++j)
            {
              for (int i = 1; i < n; ++i)
              {
                if ((j <= m / 2) && (i <= n / 2))
                  continue;

                hs[i, j] = -r[i, j];
              }
            }
          }
          else
          {
            double betta1 = 0;
            double betta2 = 0;
            for (int j = 1; j < m; ++j)
            {
              for (int i = 1; i < n; ++i)
              {
                if ((j <= m / 2) && (i <= n / 2))
                  continue;

                betta1 += ((hs[i + 1, j] - 2 * hs[i, j] + hs[i - 1, j]) * h2 + (hs[i, j + 1] - 2 * hs[i, j] + hs[i, j - 1]) * k2) * r[i, j];
                betta2 += ((hs[i + 1, j] - 2 * hs[i, j] + hs[i - 1, j]) * h2 + (hs[i, j + 1] - 2 * hs[i, j] + hs[i, j - 1]) * k2) * hs[i, j];
              }
            }
            double betta = betta1 / betta2;

            for (int j = 0; j < m + 1; ++j)
            {
              for (int i = 0; i < n + 1; ++i)
              {
                hs[i, j] = -r[i, j] + betta * hs[i, j];
              }
            }
          }

          double alpha1 = 0;
          double alpha2 = 0;
          for (int j = 1; j < m; ++j)
          {
            for (int i = 1; i < n; ++i)
            {
              if ((j <= m / 2) && (i <= n / 2))
                continue;

              alpha1 += r[i, j] * hs[i, j];
              alpha2 += ((hs[i + 1, j] - 2 * hs[i, j] + hs[i - 1, j]) * h2 + (hs[i, j + 1] - 2 * hs[i, j] + hs[i, j - 1]) * k2) * hs[i, j];
            }
          }

          double alpha = -alpha1 / alpha2;
          if (Math.Abs(alpha) > 1e-12)
          {
            for (int j = 1; j < m; ++j)
            {
              for (int i = 1; i < n; ++i)
              {
                if ((j <= m / 2) && (i <= n / 2))
                  continue;

                vOld = v[i, j];
                vNew = v[i, j] + alpha * hs[i, j];

                curEps = Math.Abs(vOld - vNew);
                if (curEps > currentEps)
                {
                  currentEps = curEps;
                }
                v[i, j] = vNew;

              }
            }
            numStep++;
          }

        } while (currentEps > epsMin && numStep <= maxStep);
        numStep--;
      }
  }
}
